Drilling Data Analysis & ROP Machine Learning Modeling Report

---

1. Executive Summary

This report documents the machine learning workflow for predicting Rate of Penetration (ROP) and analyzing real-time drilling dynamics using logging data from Well 58-32 (Pason log dataset).

The primary objective is to:

• Evaluate exploratory features
• Analyze dynamic interactions across drilling variables
• Train multi-model machine learning algorithms capable of adapting to downhole formation shifts

---

2. Key Results & Model Performance Summary

Core Takeaway: The machine learning modeling workflow successfully predicts Rate of Penetration (ROP) using key surface drilling parameters. Below is a snapshot of the primary findings and performance metrics across evaluated models.

Performance Overview (Holdout Test Block)

Model Architecture	(R^2) Score	MAE (ft/hr)	RMSE (ft/hr)	Status / Strategy
Tuned Gradient Boosting	0.7219	1.53	2.55	Primary Model (Optimized via TimeSeriesSplit CV)

3. Feature Discovery & Data Cleaning

The raw surface logging dataset was ingested and sanitized by stripping whitespace and isolating numeric channels.

• Dataset Source: Well_58-32_processed_pason_log.csv
• Total Identified Numeric Channels: 27 channels across primary operational parameters, mud properties, and gas monitoring.
• Unit Redundancy: Parameters include dual-unit entries (Imperial/Metric), which require feature pruning during pre-processing to eliminate collinearity.

Isolated Logging Parameters

Category	Parameters Included
Geometry & Penetration	Depth (ft/m), ROP (1 ft / 1 m)
Mechanical Drilling	Weight on Bit (k-lbs/kg), Hookload (k-lbs/kg), Rotary Speed (rpm), Surface Torque (psi/kPa)
Hydraulics & Mud	Pump Pressure (psi/kPa), Flow In (gal/min & L/min), Flow Out (%), WH Pressure (psi/kPa), Pit Total (bbls/m³)
Thermal Dynamics	Temp In (°F/°C), Temp Out (°F/°C)
Safety & Environmental	H2S Floor, H2S Cellar, H2S Pits

---

4. Exploratory Data Analysis & Strategy

We begin with exploratory data analysis to select optimal features and targets, while determining which machine learning tasks are feasible based on the target variable.

To kick off the analysis, we generate pairwise plots to evaluate the interaction of individual features with each other

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

file_path = r"C:\Users\USER\Downloads\Well_58-32_processed_pason_log.csv"
df = pd.read_csv(file_path)
# 1. Select only the key interacting drilling features
key_features = [
    'Depth (ft)', 
    'Weight on Bit (k-lbs)', 
    'Rotary Speed (rpm)', 
    'Surface Torque (psi)', 
    'Pump Pressure (psi)', 
    'ROP'  # Target variable
]

# Filter dataset for available key features (handles variations in column naming)
selected_cols = [col for col in key_features if col in df.columns]
df_subset = df[selected_cols].dropna()

# 2. Configure Seaborn styling for inline report rendering
sns.set_theme(style="ticks", palette="muted")

# 3. Generate focused Pairwise Interaction Plot
g = sns.pairplot(
    df_subset, 
    diag_kind="kde",           # Smooth kernel density estimate on the diagonal
    plot_kws={'alpha': 0.5, 's': 15, 'edgecolor': 'none'}, # Semi-transparent scatter points
    height=2.0,                 # Compact height per subplot to keep image size manageable
    corner=True                 # Displays only lower triangle to eliminate duplicate plots
)

g.fig.suptitle("Key Drilling Parameter Interaction Matrix", y=1.02, fontsize=14, fontweight='bold')

# 4. Display immediately inline in the notebook
plt.tight_layout()
plt.show()

2.2 Pairwise & Interaction Analysis

To examine feature-to-feature relationships, a pairwise analysis matrix covering 351 unique parameter combinations was generated (Drilling_Pairwise_Plots.pdf).

• Key Observations: Multi-modal clusters in Weight on Bit (WOB) vs. ROP and Torque vs. RPM indicate distinct lithology transitions downhole.
• Feature Pruning Recommendation: Dual metric/imperial channels (e.g., Depth(ft) vs Depth(m)) must be reduced to single features to avoid perfect multicollinearity ($r = 1.00$).

---

3. Workflow Implementation Strategy

  [ Data Ingestion & Sanitization ]
                  │
                  ▼
  [ Exploratory Pairwise Analysis ]
         │                 │
         ▼                 ▼
[ Time-Series Dynamics ]  [ Correlation Matrix ]
         │                 │
         └────────┬────────┘
                  │
                  ▼
   [ Feature Engineering & Scaling ]
                  │
                  ▼
 [ Multi-Model Training & Evaluation ]
(SVR, Random Forest, KNN, PolyReg,
 Gradient Boost, MLP Neural Net)
                  │
                  ▼
[ Depth Track Overlay & Downhole Validation ]

---

We do a time-series plot to see how the individual features interact with time.

import pandas as pd
import matplotlib.pyplot as plt

# --- 1. AUTOMATIC DATAFRAME FALLBACK ---
if 'df' not in locals():
    possible_dfs = [v for k, v in list(locals().items()) if isinstance(v, pd.DataFrame)]
    df = possible_dfs[-1] if possible_dfs else None

# --- 2. FLEXIBLE COLUMN MATCHING ---
# Map common column name variations to clean display labels
channel_mapping = {
    'ROP(1 ft)': 'ROP', 'ROP (ft/hr)': 'ROP', 'ROP': 'ROP',
    'WOB(k-lbs)': 'WOB (k-lbs)', 'Weight on Bit (k-lbs)': 'WOB (k-lbs)', 'WOB': 'WOB (k-lbs)',
    'Rotary Speed(rpm)': 'RPM', 'Rotary Speed (rpm)': 'RPM', 'RPM': 'RPM',
    'Pump Press(psi)': 'Pump Press (psi)', 'Pump Pressure (psi)': 'Pump Press (psi)', 'Pump Press': 'Pump Press (psi)'
}

# Find matching columns in your DataFrame
found_cols = [col for col in df.columns if col in channel_mapping]

if not found_cols:
    # Fallback to numeric columns if named columns aren't matched
    found_cols = list(df.select_dtypes(include='number').columns[:4])

# Remove duplicates while keeping mapping
matched_labels = {col: channel_mapping.get(col, col) for col in found_cols}
n_plots = len(matched_labels)

# --- 3. CREATE SUBPLOTS (HANDLES 1 OR MULTIPLE PLOTS) ---
fig, axes = plt.subplots(n_plots, 1, figsize=(10, max(2.5, 2.2 * n_plots)), sharex=True)

# Ensure axes is always iterable (fixes the TypeError when n_plots == 1)
if n_plots == 1:
    axes = [axes]

# Determine X-axis (Depth channel or Index)
depth_col = next((col for col in df.columns if 'depth' in col.lower()), None)
x_axis = df[depth_col] if depth_col else df.index
x_label = depth_col if depth_col else 'Sequential Log Index'

# --- 4. PLOT EACH CHANNEL ---
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd']

for i, (col, label) in enumerate(matched_labels.items()):
    axes[i].plot(x_axis, df[col], color=colors[i % len(colors)], linewidth=1, alpha=0.85)
    axes[i].set_ylabel(label, fontsize=10, fontweight='bold')
    axes[i].grid(True, linestyle=':', alpha=0.6)

axes[-1].set_xlabel(x_label, fontsize=11, fontweight='bold')
fig.suptitle("Real-Time Surface Drilling Log Dynamics", y=0.99, fontsize=13, fontweight='bold')

plt.tight_layout()
plt.show()

Correlation Matrix & Multi-Variable Relationships

• Overview: Generated a correlation matrix to measure pairwise linear correlation coefficients ($r$) across all active drilling operational variables in the dataset.

import os
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

# 1. Load dataset
file_path = r"C:\Users\USER\Downloads\Well_58-32_processed_pason_log.csv"
df = pd.read_csv(file_path) # or .csv

if file_path.endswith(".xlsx") or file_path.endswith(".xls"):
  df = pd.read_excel(file_path)
else:
  df = pd.read_csv(file_path)

# Strip spaces from header names
df.columns = df.columns.astype(str).str.strip()

# Print actual available columns in PyCharm console to verify
print("Actual columns found in your Excel file:")
print(df.columns.tolist())

# 2. Flexible keyword search: finds columns containing these terms
keywords = ["Depth", "ROP", "weight", "Temp Out", "Temp In", "Pit", "Pump"]

matched_cols = []
for kw in keywords:
  for col in df.columns:
    if kw.lower() in col.lower() and col not in matched_cols:
      matched_cols.append(col)

print(f"\nSuccessfully matched {len(matched_cols)} columns:")
print(matched_cols)

# 3. Clean numeric data
df_subset = df[matched_cols].apply(pd.to_numeric, errors="coerce")

# Drop any column that ended up completely blank
df_subset = df_subset.dropna(how="all", axis=1)

# 4. Calculate Correlation Matrix
corr_matrix = df_subset.corr(method="pearson")

# 5. Plot Heatmap
plt.figure(figsize=(10, 8))

sns.heatmap(
    corr_matrix,
    annot=True,
    fmt=".2f",
    cmap="vlag",
    vmax=1,
    vmin=-1,
    linewidths=0.8,
    cbar_kws={"label": "Correlation Coefficient"},
)

plt.title(
    "Well 58-32 Log - Correlation Matrix (Matched Variables)",
    fontsize=13,
    fontweight="bold",
    pad=12,
)
plt.xticks(rotation=45, ha="right", fontsize=9)
plt.yticks(fontsize=9)
plt.tight_layout()

# Save image
output_file = "Well_58-32_Correlation_Matrix.png"
plt.savefig(output_file, dpi=300)
plt.show()

print(f"\nFull heatmap saved as '{output_file}'.")

Actual columns found in your Excel file:
['Depth(ft)', 'Depth(m)', 'ROP(1 ft)', 'ROP(1 m)', 'weight on bit (k-lbs)', 'weight on bit (kg)', 'Temp Out( degF)', 'Temp Out( degC)', 'Temp In(degF)', 'Temp In(degC)', 'Pit Total (bbls)', 'Pit Total (m3)', 'Pump Press (psi)', 'Pump Press (KPa)', 'Hookload (k-lbs)', 'Hookload (kg)', 'Surface Torque (psi)', 'Surface Torque (KPa)', 'Rotary Speed (rpm)', 'Flow In (gal/min)', 'Flow In(liters/min)', 'Flow Out %', 'WH Pressure (psi)', 'WH Pressure (KPa)', 'H2S Floor', 'H2S Cellar', 'H2S Pits']

Successfully matched 15 columns:
['Depth(ft)', 'Depth(m)', 'ROP(1 ft)', 'ROP(1 m)', 'weight on bit (k-lbs)', 'weight on bit (kg)', 'Temp Out( degF)', 'Temp Out( degC)', 'Temp In(degF)', 'Temp In(degC)', 'Pit Total (bbls)', 'Pit Total (m3)', 'H2S Pits', 'Pump Press (psi)', 'Pump Press (KPa)']

Full heatmap saved as 'Well_58-32_Correlation_Matrix.png'.

Phase 1: Data Ingestion & Preprocessing Workflow

---

Overview & Next Steps

Having completed the exploratory data analysis and correlation matrix evaluation across the full dataset, we now begin the machine learning model development pipelinespan_1span_1.

• Objective: Ingest, sanitize, and clean the dataset to construct a clean feature matrix for downstream model trainingspan_2span_2.

import pandas as pd
import numpy as np

pd.set_option('display.max_columns', None)
pd.set_option('display.width', 1000)

print("=" * 70)
print("PHASE 1: DATA LOADING, CLEANING & DRILLING FILTERING")
print("=" * 70)

# 1. LOAD DATA (Update path as needed)
file_path = r"C:\Users\USER\Downloads\Well_58-32_processed_pason_log.csv"
print("\n[STEP 1] Loading Drilling Log File...")
df = pd.read_csv(file_path)
print(f"-> Raw File Loaded: {df.shape[0]} rows across {df.shape[1]} columns.")

# 2. CHRONOLOGICAL SORTING BY DEPTH
print("\n" + "-" * 70)
print("[STEP 2] Sorting Chronologically / Depth-wise...")
depth_col = 'Depth(ft)' if 'Depth(ft)' in df.columns else 'Depth'
df = df.sort_values(depth_col).reset_index(drop=True)

# 3. FILTER ACTIVE DRILLING INTERVALS
print("\n" + "-" * 70)
print("[STEP 3] Filtering Active Drilling Days/Intervals...")
target_col = 'ROP(1 ft)' if 'ROP(1 ft)' in df.columns else 'ROP'

# Retain only active drilling intervals (filtering out connections/tripping)
active_mask = (
    (df[target_col] > 0) & 
    (df['weight on bit (k-lbs)'] > 0) & 
    (df['Rotary Speed (rpm)'] > 0)
)
active_df = df[active_mask].copy()

print(f"-> Total Raw Log Records: {len(df)}")
print(f"-> Active Drilling Intervals Preserved: {len(active_df)}")
print(f"-> Non-Drilling / Connection / Tripping Noise Filtered: {len(df) - len(active_df)}")

# 4. SELECT KEY OPERATIONAL FEATURES (Deduplicating Unit Pairs)
print("\n" + "-" * 70)
print("[STEP 4] Selecting Target & Single-Unit Features...")

candidate_features = [
    depth_col,
    'weight on bit (k-lbs)',
    'Rotary Speed (rpm)',
    'Pump Press (psi)',
    'Flow In (gal/min)',
    'Surface Torque (psi)',
    'Hookload (k-lbs)',
    'Temp In(degF)',
    'Temp Out( degF)',
    target_col
]

selected_cols = [c for c in candidate_features if c in active_df.columns]
ml_subset = active_df[selected_cols].copy()

print(f"-> Columns Preserved ({len(selected_cols)} total):\n   {selected_cols}")

# 5. EXPORT CLEANED FILE FOR PHASE 2
output_filename = "drilling_phase1_cleaned.csv"
ml_subset.to_csv(output_filename, index=False)
print(f"\n-> Phase 1 Output Saved to '{output_filename}'")

print("\n" + "=" * 70)
print("PHASE 1 COMPLETE: RECOVERED ACTIVE DRILLING HISTORY")
print("=" * 70)

======================================================================
PHASE 1: DATA LOADING, CLEANING & DRILLING FILTERING
======================================================================

[STEP 1] Loading Drilling Log File...
-> Raw File Loaded: 7311 rows across 27 columns.

----------------------------------------------------------------------
[STEP 2] Sorting Chronologically / Depth-wise...

----------------------------------------------------------------------
[STEP 3] Filtering Active Drilling Days/Intervals...
-> Total Raw Log Records: 7311
-> Active Drilling Intervals Preserved: 6814
-> Non-Drilling / Connection / Tripping Noise Filtered: 497

----------------------------------------------------------------------
[STEP 4] Selecting Target & Single-Unit Features...
-> Columns Preserved (10 total):
   ['Depth(ft)', 'weight on bit (k-lbs)', 'Rotary Speed (rpm)', 'Pump Press (psi)', 'Flow In (gal/min)', 'Surface Torque (psi)', 'Hookload (k-lbs)', 'Temp In(degF)', 'Temp Out( degF)', 'ROP(1 ft)']

-> Phase 1 Output Saved to 'drilling_phase1_cleaned.csv'

======================================================================
PHASE 1 COMPLETE: RECOVERED ACTIVE DRILLING HISTORY
======================================================================

Phase 2: Feature Engineering & Temporal Encoding

---

Phase 1 Completion Summary

• Saved Dataset: drilling_phase1_cleaned.csv
• Preserved Channels (10 Total): Depth(ft), Weight on Bit (k-lbs), Rotary Speed (rpm), Pump Press (psi), Flow In (gal/min), Surface Torque (psi), Temp In (degF), Temp Out (degF), ROP(1 ft)

---

Phase 2 Overview

With Phase 1 complete and active drilling history recovered, we move into feature engineering, temporal encoding, and integrating domain drilling physics to prepare the dataset for predictive modeling.

import pandas as pd
import numpy as np

pd.set_option('display.max_columns', None)
pd.set_option('display.width', 1000)

print("=" * 70)
print("PHASE 2: FEATURE ENGINEERING & DRILLING PHYSICS")
print("=" * 70)

# 0. LOAD PHASE 1 CLEANED DATA
file_input = "drilling_phase1_cleaned.csv"
df_p2 = pd.read_csv(file_input)

depth_col = 'Depth(ft)' if 'Depth(ft)' in df_p2.columns else 'Depth'
target_col = 'ROP(1 ft)' if 'ROP(1 ft)' in df_p2.columns else 'ROP'

df_p2 = df_p2.sort_values(depth_col).reset_index(drop=True)
print(f"\n[STEP 0] Base Data Loaded: {len(df_p2)} rows across {df_p2.shape[1]} columns.")

# 1. SENSOR IMPUTATION
print("\n" + "-" * 70)
print("[STEP 1] Imputing Missing Sensor Readings...")
sensor_cols = [c for c in df_p2.columns if c not in [depth_col, target_col]]
df_p2[sensor_cols] = df_p2[sensor_cols].ffill().bfill()
print("-> Successfully forward/backward filled missing operational sensors.")

# 2. LAGS & ROLLING TRENDS
print("\n" + "-" * 70)
print("[STEP 2] Calculating Historical Lags & Rolling Averages...")
for col in sensor_cols:
    df_p2[f'{col}_lag1'] = df_p2[col].shift(1)
    df_p2[f'{col}_lag2'] = df_p2[col].shift(2)
    df_p2[f'{col}_roll3_mean'] = df_p2[col].shift(1).rolling(window=3, min_periods=1).mean()
    df_p2[f'{col}_roll7_mean'] = df_p2[col].shift(1).rolling(window=7, min_periods=1).mean()

df_p2['TARGET_lag1'] = df_p2[target_col].shift(1)
df_p2['TARGET_lag2'] = df_p2[target_col].shift(2)

# 3. PHYSICAL INTERACTION TERMS
print("\n" + "-" * 70)
print("[STEP 3] Engineering Drilling Physics Interactions...")

wob_col = 'weight on bit (k-lbs)'
rpm_col = 'Rotary Speed (rpm)'
flow_col = 'Flow In (gal/min)'
press_col = 'Pump Press (psi)'

if wob_col in df_p2.columns and rpm_col in df_p2.columns:
    df_p2['WOB_RPM_INTERACTION'] = df_p2[wob_col] * df_p2[rpm_col]

if flow_col in df_p2.columns and press_col in df_p2.columns:
    df_p2['HYDRAULIC_POWER_PROXY'] = df_p2[flow_col] * df_p2[press_col]

if 'Temp Out( degF)' in df_p2.columns and 'Temp In(degF)' in df_p2.columns:
    df_p2['DELTA_TEMP'] = df_p2['Temp Out( degF)'] - df_p2['Temp In(degF)']

# 4. REMOVE BOUNDARY LAG ROWS
print("\n" + "-" * 70)
print("[STEP 4] Removing Boundary Window NaNs...")
initial_rows = len(df_p2)
df_p2_clean = df_p2.dropna(subset=['TARGET_lag2']).reset_index(drop=True)
dropped_rows = initial_rows - len(df_p2_clean)

print(f"-> Total Rows Before Cleanup: {initial_rows}")
print(f"-> Initial Window Rows Removed: {dropped_rows}")
print(f"-> Final Operational Feature Matrix: {df_p2_clean.shape[0]} rows x {df_p2_clean.shape[1]} columns")

# EXPORT
output_filename_p2 = "drilling_phase2_engineered.csv"
df_p2_clean.to_csv(output_filename_p2, index=False)
print(f"\n-> Phase 2 Output Saved to '{output_filename_p2}'")

print("\n" + "=" * 70)
print("PHASE 2 COMPLETE: READY FOR MODELING")
print("=" * 70)

======================================================================
PHASE 2: FEATURE ENGINEERING & DRILLING PHYSICS
======================================================================

[STEP 0] Base Data Loaded: 6814 rows across 10 columns.

----------------------------------------------------------------------
[STEP 1] Imputing Missing Sensor Readings...
-> Successfully forward/backward filled missing operational sensors.

----------------------------------------------------------------------
[STEP 2] Calculating Historical Lags & Rolling Averages...

----------------------------------------------------------------------
[STEP 3] Engineering Drilling Physics Interactions...

----------------------------------------------------------------------
[STEP 4] Removing Boundary Window NaNs...
-> Total Rows Before Cleanup: 6814
-> Initial Window Rows Removed: 2
-> Final Operational Feature Matrix: 6812 rows x 47 columns

-> Phase 2 Output Saved to 'drilling_phase2_engineered.csv'

======================================================================
PHASE 2 COMPLETE: READY FOR MODELING
======================================================================

4. Machine Learning Modeling & Architecture

---

4.1 Candidate Models Evaluated

Six distinct algorithmic architectures were configured to predict continuous ROP values:

1. Support Vector Machine (SVM): Kernel-based regression for non-linear decision boundaries.
2. Random Forest: Ensemble decision trees capturing non-linear feature interaction with resilience to outliers.
3. K-Nearest Neighbors (KNN): Non-parametric spatial distance estimation across drilling regimes.
4. Polynomial Regression: Baseline multi-variable linear model with polynomial interaction terms.
5. Gradient Boosted Decision Trees (GBDT): Sequential boosting targeting residual error minimization.
6. Multilayer Perceptron (MLP): Deep neural network architecture designed for complex downhole interactions.

---

4.2 Target & Feature Selection

• Target Variable: Rate of Penetration (ROP)
• Primary Feature Predictors: Depth, Weight on Bit (WOB), Rotary Speed (RPM), Surface Torque, Flow In, Pump Pressure

---

We evaluate Support Vector Machine, Random Forest, K-Nearest Neighbors, Polynomial Regression, Gradient Boosted Decision Trees, and Multilayer Perceptron (Neural Net) model performance.

import pandas as pd
import numpy as np

from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.svm import SVR
from sklearn.neighbors import KNeighborsRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

pd.set_option('display.max_columns', None)
pd.set_option('display.width', 1000)

print("=" * 70)
print("PHASE 3: MULTI-MODEL ML BENCHMARK (ROP PREDICTION)")
print("=" * 70)

# 0. LOAD PHASE 2 ENGINEERED DATA
file_input = "drilling_phase2_engineered.csv"
df_p3 = pd.read_csv(file_input)

depth_col = 'Depth(ft)' if 'Depth(ft)' in df_p3.columns else 'Depth'
target_col = 'ROP(1 ft)' if 'ROP(1 ft)' in df_p3.columns else 'ROP'

# Do NOT ignore depth; only ignore the target column during feature selection
ignore_cols = [target_col]
feature_cols = [col for col in df_p3.columns if col not in ignore_cols]

# CLEANING STEP: Handle Inf/NaNs
df_p3[feature_cols] = df_p3[feature_cols].replace([np.inf, -np.inf], np.nan)
df_p3 = df_p3.dropna(subset=feature_cols + [target_col]).reset_index(drop=True)

print(f"\n[STEP 0] Cleaned Feature Matrix Loaded: {len(df_p3)} active rows.")
print(f"-> Predictor Features: {len(feature_cols)}")

from sklearn.model_selection import train_test_split

from sklearn.utils import shuffle

# 1. BLOCK SPLIT FOLLOWED BY TRAIN-SET SHUFFLING
print("\n" + "-" * 70)
print("[STEP 1] Performing Block Depth Split then Shuffling Train Set...")

# Step 1A: Sort chronologically by depth to establish the continuous profile
df_p3 = df_p3.sort_values(depth_col).reset_index(drop=True)

# Step 1B: Calculate 80/20 depth cutoff index
split_idx = int(len(df_p3) * 0.80)

# Step 1C: Slice into shallow (train) and deep (test) continuous blocks
train_block = df_p3.iloc[:split_idx]
test_block  = df_p3.iloc[split_idx:]

# Step 1D: Shuffle ONLY the training block internally
train_block_shuffled = shuffle(train_block, random_state=42).reset_index(drop=True)

# Step 1E: Assign features and target
X_train = train_block_shuffled[feature_cols]
y_train = train_block_shuffled[target_col]

# Test set remains sequential for downhole evaluation
X_test = test_block[feature_cols]
y_test = test_block[target_col]

print(f"-> Training Set (Shuffled):   {len(X_train)} samples")
print(f"-> Testing Set (Sequential): {len(X_test)} samples")
# 2. FEATURE SCALING
scaler = StandardScaler()
X_train_scaled = pd.DataFrame(scaler.fit_transform(X_train), columns=feature_cols)
X_test_scaled = pd.DataFrame(scaler.transform(X_test), columns=feature_cols)

# 3. BENCHMARK ALL 6 REQUESTED ML MODELS
print("\n" + "-" * 70)
print("[STEP 2] Training Benchmark Machine Learning Models...")
random_state = 42
models = {
    "Random Forest": (RandomForestRegressor(n_estimators=100, random_state=42), False),
    "Gradient Boosted Trees": (GradientBoostingRegressor(n_estimators=100, learning_rate=0.05, random_state=42), False),
    "Support Vector Machine (SVR)": (SVR(kernel='rbf', C=100, epsilon=0.1), True),
    "Nearest Neighbors (k-NN)": (KNeighborsRegressor(n_neighbors=5), True),
    "Polynomial Ridge Regression": (make_pipeline(PolynomialFeatures(degree=2, interaction_only=True, include_bias=False), Ridge(alpha=100.0)), True),
    "Multilayer Perceptron (Neural Net)": (
        MLPRegressor(
            hidden_layer_sizes=(64, 32),
            activation='relu',
            solver='adam',
            max_iter=1000,
            early_stopping=True,
            random_state= random_state
        ), 
        True
    )
}

results = []

for name, (model, requires_scaling) in models.items():
    X_tr = X_train_scaled if requires_scaling else X_train
    X_te = X_test_scaled if requires_scaling else X_test
    
    print(f"-> Fitting model: {name}...")
    model.fit(X_tr, y_train)
    y_pred = model.predict(X_te)
    
    r2 = r2_score(y_test, y_pred)
    rmse = np.sqrt(mean_squared_error(y_test, y_pred))
    mae = mean_absolute_error(y_test, y_pred)
    
    results.append({
        "Model": name,
        "R2 Score": round(r2, 4),
        "RMSE": round(rmse, 2),
        "MAE": round(mae, 2)
    })

# 4. SUMMARY EVALUATION TABLE
print("\n" + "-" * 70)
print("[STEP 3] MODEL PERFORMANCE COMPARISON SUMMARY")
print("-" * 70)

summary_df = pd.DataFrame(results).sort_values(by="R2 Score", ascending=False).reset_index(drop=True)
print(summary_df.to_string(index=False))

# 5. FEATURE IMPORTANCES
print("\n" + "-" * 70)
print("[STEP 4] Top 10 Feature Importances (Random Forest)")
print("-" * 70)

rf_model = models["Random Forest"][0]
importances = pd.DataFrame({
    'Feature': feature_cols,
    'Importance': rf_model.feature_importances_
}).sort_values(by='Importance', ascending=False).reset_index(drop=True)

print(importances.head(10).to_string(index=False))

print("\n" + "=" * 70)
print("PHASE 3 COMPLETE: MULTI-MODEL BENCHMARK FINISHED")
print("=" * 70)

======================================================================
PHASE 3: MULTI-MODEL ML BENCHMARK (ROP PREDICTION)
======================================================================

[STEP 0] Cleaned Feature Matrix Loaded: 6812 active rows.
-> Predictor Features: 46

----------------------------------------------------------------------
[STEP 1] Performing Block Depth Split then Shuffling Train Set...
-> Training Set (Shuffled):   5449 samples
-> Testing Set (Sequential): 1363 samples

----------------------------------------------------------------------
[STEP 2] Training Benchmark Machine Learning Models...
-> Fitting model: Random Forest...
-> Fitting model: Gradient Boosted Trees...
-> Fitting model: Support Vector Machine (SVR)...
-> Fitting model: Nearest Neighbors (k-NN)...
-> Fitting model: Polynomial Ridge Regression...
-> Fitting model: Multilayer Perceptron (Neural Net)...

----------------------------------------------------------------------
[STEP 3] MODEL PERFORMANCE COMPARISON SUMMARY
----------------------------------------------------------------------
                             Model  R2 Score  RMSE  MAE
            Gradient Boosted Trees    0.5333  3.31 2.08
                     Random Forest    0.4140  3.71 2.59
Multilayer Perceptron (Neural Net)    0.0905  4.62 3.45
          Nearest Neighbors (k-NN)   -0.7989  6.49 4.99
      Support Vector Machine (SVR)   -1.8208  8.13 5.54
       Polynomial Ridge Regression   -3.0096  9.70 6.31

----------------------------------------------------------------------
[STEP 4] Top 10 Feature Importances (Random Forest)
----------------------------------------------------------------------
                      Feature  Importance
                    Depth(ft)    0.646782
                  TARGET_lag1    0.110127
            Flow In (gal/min)    0.099153
Rotary Speed (rpm)_roll3_mean    0.034847
  Hookload (k-lbs)_roll3_mean    0.009327
                  TARGET_lag2    0.009055
       Flow In (gal/min)_lag2    0.007222
 Flow In (gal/min)_roll7_mean    0.006031
             Hookload (k-lbs)    0.005689
 Flow In (gal/min)_roll3_mean    0.005391

======================================================================
PHASE 3 COMPLETE: MULTI-MODEL BENCHMARK FINISHED
======================================================================

Phase 3: Multi-Model Benchmark & Spatial Leakage Analysis

---

From the above code, we used a depth cut first, then shuffled the training set, and it yielded the R2 scores above. We also tried using a fully randomized split, which yielded an R2 > 0.88 due to spatial autocorrelation. The block split scores of ~0.53 represents the true unleaked generalization capability of the models when penetrating undrilled rock formations downhole. Now, we want to remove depth as a feature for evaluating ROP because of depth being a spatial variable for evaluating ROP.

---

Key Takeaways for Presentation

• Shuffled Split & Spatial Autocorrelation: A fully randomized train/test split gave artificially high performance ($R^2 > 0.88$) because neighboring depth points are highly correlated.
• True Generalization (Block Split): Using a contiguous block split reduced the score to $R^2 \approx 0.53$, showing how the models actually perform on unseen, undrilled formations.
• Removing Depth: Depth was removed from the feature set because it acts as a spatial position marker rather than a dynamic physical force driving instantaneous rate of penetration.

import pandas as pd
import numpy as np

from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.svm import SVR
from sklearn.neighbors import KNeighborsRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
from sklearn.utils import shuffle

pd.set_option('display.max_columns', None)
pd.set_option('display.width', 1000)

print("=" * 70)
print("PHASE 3: MULTI-MODEL ML BENCHMARK (ROP PREDICTION - EXCLUDING DEPTH)")
print("=" * 70)

# 0. LOAD PHASE 2 ENGINEERED DATA
file_input = "drilling_phase2_engineered.csv"
df_p3 = pd.read_csv(file_input)

depth_col = 'Depth(ft)' if 'Depth(ft)' in df_p3.columns else 'Depth'
target_col = 'ROP(1 ft)' if 'ROP(1 ft)' in df_p3.columns else 'ROP'

# EXCLUDE DEPTH FROM PREDICTOR FEATURES
# (Depth is placed in ignore_cols so it isn't trained on, but remains available for sorting/splitting)
ignore_cols = [depth_col, target_col]
feature_cols = [col for col in df_p3.columns if col not in ignore_cols]

# CLEANING STEP: Handle Inf/NaNs
df_p3[feature_cols] = df_p3[feature_cols].replace([np.inf, -np.inf], np.nan)
df_p3 = df_p3.dropna(subset=feature_cols + [target_col]).reset_index(drop=True)

print(f"\n[STEP 0] Cleaned Feature Matrix Loaded: {len(df_p3)} active rows.")
print(f"-> Predictor Features (Excluding Depth): {len(feature_cols)}")

# 1. BLOCK DEPTH SPLIT FOLLOWED BY TRAIN-SET SHUFFLING
print("\n" + "-" * 70)
print("[STEP 1] Performing Block Depth Split then Shuffling Train Set...")

# Step 1A: Sort chronologically by depth to establish continuous wellbore profile
df_p3 = df_p3.sort_values(depth_col).reset_index(drop=True)

# Step 1B: Calculate 80/20 depth cutoff index
split_idx = int(len(df_p3) * 0.80)

# Step 1C: Slice into shallow (train) and deep (test) continuous blocks
train_block = df_p3.iloc[:split_idx]
test_block  = df_p3.iloc[split_idx:]

# Step 1D: Shuffle ONLY the training block internally
train_block_shuffled = shuffle(train_block, random_state=42).reset_index(drop=True)

# Step 1E: Assign features and target
X_train = train_block_shuffled[feature_cols]
y_train = train_block_shuffled[target_col]

# Test set remains sequential for true downhole evaluation
X_test = test_block[feature_cols]
y_test = test_block[target_col]

print(f"-> Training Set (Shuffled):   {len(X_train)} samples")
print(f"-> Testing Set (Sequential): {len(X_test)} samples")

# 2. FEATURE SCALING
scaler = StandardScaler()
X_train_scaled = pd.DataFrame(scaler.fit_transform(X_train), columns=feature_cols)
X_test_scaled = pd.DataFrame(scaler.transform(X_test), columns=feature_cols)

# 3. BENCHMARK ALL 6 REQUESTED ML MODELS
print("\n" + "-" * 70)
print("[STEP 2] Training Benchmark Machine Learning Models...")

models = {
    "Random Forest": (RandomForestRegressor(n_estimators=100, random_state=42), False),
    "Gradient Boosted Trees": (GradientBoostingRegressor(n_estimators=100, learning_rate=0.05, random_state=42), False),
    "Support Vector Machine (SVR)": (SVR(kernel='rbf', C=100, epsilon=0.1), True),
    "Nearest Neighbors (k-NN)": (KNeighborsRegressor(n_neighbors=5), True),
    "Polynomial Ridge Regression": (make_pipeline(PolynomialFeatures(degree=2, interaction_only=True, include_bias=False), Ridge(alpha=100.0)), True),
    "Multilayer Perceptron (Neural Net)": (
        MLPRegressor(
            hidden_layer_sizes=(64, 32),
            activation='relu',
            solver='adam',
            max_iter=1000,
            early_stopping=True,
            random_state= random_state
        ), 
        True
    )
}

results = []

for name, (model, requires_scaling) in models.items():
    X_tr = X_train_scaled if requires_scaling else X_train
    X_te = X_test_scaled if requires_scaling else X_test
    
    print(f"-> Fitting model: {name}...")
    model.fit(X_tr, y_train)
    y_pred = model.predict(X_te)
    
    r2 = r2_score(y_test, y_pred)
    rmse = np.sqrt(mean_squared_error(y_test, y_pred))
    mae = mean_absolute_error(y_test, y_pred)
    
    results.append({
        "Model": name,
        "R2 Score": round(r2, 4),
        "RMSE": round(rmse, 2),
        "MAE": round(mae, 2)
    })

# 4. SUMMARY EVALUATION TABLE
print("\n" + "-" * 70)
print("[STEP 3] MODEL PERFORMANCE COMPARISON SUMMARY")
print("-" * 70)

summary_df = pd.DataFrame(results).sort_values(by="R2 Score", ascending=False).reset_index(drop=True)
print(summary_df.to_string(index=False))

# 5. FEATURE IMPORTANCES
print("\n" + "-" * 70)
print("[STEP 4] Top 10 Feature Importances (Random Forest)")
print("-" * 70)

rf_model = models["Random Forest"][0]
importances = pd.DataFrame({
    'Feature': feature_cols,
    'Importance': rf_model.feature_importances_
}).sort_values(by='Importance', ascending=False).reset_index(drop=True)

print(importances.head(10).to_string(index=False))

print("\n" + "=" * 70)
print("PHASE 3 COMPLETE: MULTI-MODEL BENCHMARK FINISHED")
print("=" * 70)

======================================================================
PHASE 3: MULTI-MODEL ML BENCHMARK (ROP PREDICTION - EXCLUDING DEPTH)
======================================================================

[STEP 0] Cleaned Feature Matrix Loaded: 6812 active rows.
-> Predictor Features (Excluding Depth): 45

----------------------------------------------------------------------
[STEP 1] Performing Block Depth Split then Shuffling Train Set...
-> Training Set (Shuffled):   5449 samples
-> Testing Set (Sequential): 1363 samples

----------------------------------------------------------------------
[STEP 2] Training Benchmark Machine Learning Models...
-> Fitting model: Random Forest...
-> Fitting model: Gradient Boosted Trees...
-> Fitting model: Support Vector Machine (SVR)...
-> Fitting model: Nearest Neighbors (k-NN)...
-> Fitting model: Polynomial Ridge Regression...
-> Fitting model: Multilayer Perceptron (Neural Net)...

----------------------------------------------------------------------
[STEP 3] MODEL PERFORMANCE COMPARISON SUMMARY
----------------------------------------------------------------------
                             Model  R2 Score  RMSE  MAE
Multilayer Perceptron (Neural Net)   -0.0324  4.92 3.85
          Nearest Neighbors (k-NN)   -0.8310  6.55 5.05
            Gradient Boosted Trees   -1.3525  7.43 2.88
      Support Vector Machine (SVR)   -1.4128  7.52 4.78
                     Random Forest   -2.4223  8.96 3.49
       Polynomial Ridge Regression   -5.6911 12.53 8.64

----------------------------------------------------------------------
[STEP 4] Top 10 Feature Importances (Random Forest)
----------------------------------------------------------------------
                      Feature  Importance
  Pump Press (psi)_roll7_mean    0.574184
            Flow In (gal/min)    0.120513
                  TARGET_lag1    0.115452
Rotary Speed (rpm)_roll3_mean    0.050292
       Flow In (gal/min)_lag2    0.015635
                  TARGET_lag2    0.011685
 Flow In (gal/min)_roll7_mean    0.011233
  Hookload (k-lbs)_roll3_mean    0.008932
             Pump Press (psi)    0.007253
  Pump Press (psi)_roll3_mean    0.006902

======================================================================
PHASE 3 COMPLETE: MULTI-MODEL BENCHMARK FINISHED
======================================================================

Block Split Dynamics & Mechanical Specific Energy (MSE) Introduction

---

Original Process Commentary

When we performed a Block Split (train on the top 80% of the well, test on the bottom 20%), the test set is in a deeper, more compacted formation that the model has never seenspan_2span_2.

Depth acts as the key geological anchor: Subsurface compaction, pore pressure, and rock hardness change non-linearly downholespan_3span_3. Depth serves as the primary proxy for those unseen rock propertiesspan_4span_4.

Without Depth, models cannot scale to deeper formations: As seen in our Feature Importance list, Pump Press (psi)_roll17_mean became the primary driver (57.4% importance)span_5span_5. However, surface parameters like pump pressure and flow rate alone do not reflect the increased rock strength at greater depths, causing the models to systematically over-predict or under-predict ROP in the test blockspan_6span_6.

To make the model better, We tried to calculate Mechanical Specific Energy (MSE), which combines: WOB, RPM, Torque and ROP in a single parameter reflecting Rock drillability independent of Depth. with the codes belowspan_7span_7

---

Key Takeaways:

• Block Split Impact (80/20): Training on the upper 80% and testing on the lower 20% forces the model to evaluate a deeper, more compacted formation it has never encounteredspan_8span_8.
• Depth as a Geological Proxy: Subsurface compaction, pore pressure, and rock hardness change non-linearly downhole, making depth the main proxy for these unseen rock propertiesspan_9span_9.
• Feature Shift without Depth: Removing depth shifted primary importance to Pump Press (psi)_roll17_mean (57.4%), but surface parameters alone cannot capture rock strength changes at deeper levelsspan_10span_10.
• MSE Integration: Calculating Mechanical Specific Energy (MSE) combines WOB, RPM, Torque, and ROP into a unified metric to represent rock drillability independently of depthspan_11span_11.

4 INTEGRATION OF MECHANICAL SPECIFIC ENERGY (MSE)

Core Objective: MSE combines operational variables (WOB, RPM, Torque, and ROP) into a single unified parameter that quantifies rock drillability and mechanical efficiency independently of depth.

The implementation and feature calculation details are executed in the code block below:

4.1 Physics-Informed Feature Engineering: Mechanical Specific Energy (MSE)

To eliminate dependency on raw Depth—which causes algorithms to overfit to specific depth intervals rather than learning true formation drillability—we engineered Mechanical Specific Energy (MSE).

---

Mathematical Formulation

$$\text{MSE} = \frac{F_B}{A_B} + \frac{120 \pi \cdot N \cdot T}{A_B \cdot \text{ROP}}$$

Where:

• $F_B$ = Weight on Bit ($\text{lbs}$)
• $A_B$ = Cross-sectional area of the drill bit ($\text{in}^2$)
• $N$ = Rotary Speed ($\text{RPM}$)
• $T$ = Surface Torque ($\text{ft-lbs}$)
• $\text{ROP}$ = Rate of Penetration ($\text{ft/hr}$)

---

Physical Proxy Mechanism

Because $\text{MSE}$ continuously incorporates applied mechanical energy ($WOB$, $RPM$, $Torque$) against output rate ($\text{ROP}$), it serves as a direct, depth-independent mathematical proxy for Unconfined Compressive Strength (UCS).

This enables algorithms like Gradient Boosted Decision Trees to detect downhole formation hardness transitions and adapt to varying rock strength purely through surface drilling physics, without explicitly relying on raw depth features.

---

Core Chain of Thought (Quick Review)

The Physics: $\text{MSE}$ combines axial force ($\frac{F_B}{A_B}$) and rotary work ($\frac{120 \pi \cdot N \cdot T}{A_B \cdot \text{ROP}}$) required to remove a unit volume of rock.

The Proxy Function: $\text{MSE}$ directly reflects rock hardness (Unconfined Compressive Strength / UCS).

The ML Advantage: Replacing Depth with $\text{MSE}$ lets tree models sense downhole lithology changes purely through surface sensor mechanics, preventing depth-interval overfitting.

4.2 Algorithm Failure Analysis & Architectural Diagnostics

While tree-based ensemble methods successfully adapted to the engineered features, several baseline models failed during regression evaluation:

• Gradient Boosted Decision Trees ((R^2 = 0.5390), (\text{MAE} = 1.77)): Decision trees partition feature space using orthogonal cuts, making them invariant to monotonic feature scaling and robust against collinear lagged features. The sequential boosting approach effectively minimized residual errors across complex non-linear intervals.
• Random Forest ((R^2 = 0.2782), (\text{MAE} = 2.24)): Averaged tree predictions provided stability, though decision trees slightly overfit to dominant rolling average features.

$k$-Nearest Neighbors ($R^2 = -0.2702$): $k\text{-NN}$ relies on Euclidean distance in a $47$-dimensional space. Due to the high dimensionality ("curse of dimensionality") and temporal noise, distance metrics lose discriminative power.Support Vector Regression ($R^2 = -4.0890$): SVR's Radial Basis Function (RBF) kernel struggled to define hyperplanes due to feature collinearity introduced by rolling window and lag variables.Polynomial Ridge Regression ($R^2 = -2986.77$): Second-degree interaction terms generated over $1,000$ cross-product terms across 47 features. Multicollinearity across lagged variables led to severe matrix ill-conditioning and massive variance explosion on unseen sequential test intervals.

import pandas as pd
import numpy as np

from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.svm import SVR
from sklearn.neighbors import KNeighborsRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
from sklearn.utils import shuffle

pd.set_option('display.max_columns', None)
pd.set_option('display.width', 1000)

print("=" * 70)
print("PHASE 3: MULTI-MODEL ML BENCHMARK (MSE FEATURE ENRICHED - NO DEPTH)")
print("=" * 70)

# 0. LOAD PHASE 2 ENGINEERED DATA
file_input = "drilling_phase2_engineered.csv"
df_p3 = pd.read_csv(file_input)

depth_col = 'Depth(ft)' if 'Depth(ft)' in df_p3.columns else 'Depth'
target_col = 'ROP(1 ft)' if 'ROP(1 ft)' in df_p3.columns else 'ROP'

# 0A. ENGINEER MECHANICAL SPECIFIC ENERGY (MSE)
# Standard 8.5-inch hole size area (sq in)
bit_diameter = 8.5
bit_area = (np.pi / 4.0) * (bit_diameter ** 2)

# Identify required columns dynamically
wob_col = [c for c in df_p3.columns if 'WOB' in c.upper() or 'WEIGHT ON BIT' in c.upper()]
rpm_col = [c for c in df_p3.columns if 'RPM' in c.upper() or 'ROTARY' in c.upper()]
torque_col = [c for c in df_p3.columns if 'TORQUE' in c.upper()]

if wob_col and rpm_col and torque_col:
    wob = df_p3[wob_col[0]]
    rpm = df_p3[rpm_col[0]]
    torque = df_p3[torque_col[0]]
    rop = df_p3[target_col].replace(0, np.nan)  # Avoid division by zero
    
    # Teale's MSE Equation (psi)
    df_p3['MSE_psi'] = (wob / bit_area) + ((120 * np.pi * rpm * torque) / (bit_area * rop))
    df_p3['MSE_psi_roll5_mean'] = df_p3['MSE_psi'].rolling(window=5, min_periods=1).mean()
    print("-> Successfully Engineered 'MSE_psi' and 'MSE_psi_roll5_mean' Features.")
else:
    print("-> Warning: WOB, RPM, or Torque columns not found. Skipping direct MSE calculation.")

# EXCLUDE DEPTH FROM PREDICTOR FEATURES
ignore_cols = [depth_col, target_col]
feature_cols = [col for col in df_p3.columns if col not in ignore_cols]

# CLEANING STEP: Handle Inf/NaNs resulting from MSE or rolling functions
df_p3[feature_cols] = df_p3[feature_cols].replace([np.inf, -np.inf], np.nan)
df_p3 = df_p3.dropna(subset=feature_cols + [target_col]).reset_index(drop=True)

print(f"\n[STEP 0] Cleaned Feature Matrix Loaded: {len(df_p3)} active rows.")
print(f"-> Predictor Features (Excluding Raw Depth): {len(feature_cols)}")

# 1. BLOCK DEPTH SPLIT FOLLOWED BY TRAIN-SET SHUFFLING
print("\n" + "-" * 70)
print("[STEP 1] Performing Block Depth Split then Shuffling Train Set...")

# Step 1A: Sort chronologically by depth to establish continuous wellbore profile
df_p3 = df_p3.sort_values(depth_col).reset_index(drop=True)

# Step 1B: Calculate 80/20 depth cutoff index
split_idx = int(len(df_p3) * 0.80)

# Step 1C: Slice into shallow (train) and deep (test) continuous blocks
train_block = df_p3.iloc[:split_idx]
test_block  = df_p3.iloc[split_idx:]

# Step 1D: Shuffle ONLY the training block internally
train_block_shuffled = shuffle(train_block, random_state=42).reset_index(drop=True)

# Step 1E: Assign features and target
X_train = train_block_shuffled[feature_cols]
y_train = train_block_shuffled[target_col]

# Test set remains sequential for true downhole evaluation
X_test = test_block[feature_cols]
y_test = test_block[target_col]

print(f"-> Training Set (Shuffled):   {len(X_train)} samples")
print(f"-> Testing Set (Sequential): {len(X_test)} samples")

# 2. FEATURE SCALING
scaler = StandardScaler()
X_train_scaled = pd.DataFrame(scaler.fit_transform(X_train), columns=feature_cols)
X_test_scaled = pd.DataFrame(scaler.transform(X_test), columns=feature_cols)

# 3. BENCHMARK ALL 6 REQUESTED ML MODELS
print("\n" + "-" * 70)
print("[STEP 2] Training Benchmark Machine Learning Models...")

models = {
    "Random Forest": (RandomForestRegressor(n_estimators=100, random_state=42), False),
    "Gradient Boosted Trees": (GradientBoostingRegressor(n_estimators=100, learning_rate=0.05, random_state=42), False),
    "Support Vector Machine (SVR)": (SVR(kernel='rbf', C=100, epsilon=0.1), True),
    "Nearest Neighbors (k-NN)": (KNeighborsRegressor(n_neighbors=5), True),
    "Polynomial Ridge Regression": (make_pipeline(PolynomialFeatures(degree=2, interaction_only=True, include_bias=False), Ridge(alpha=100.0)), True),
    "Multilayer Perceptron (Neural Net)": (
        MLPRegressor(
            hidden_layer_sizes=(64, 32),
            activation='relu',
            solver='adam',
            max_iter=1000,
            early_stopping=True,
            random_state= random_state
        ), 
        True
    )
}

results = []

for name, (model, requires_scaling) in models.items():
    X_tr = X_train_scaled if requires_scaling else X_train
    X_te = X_test_scaled if requires_scaling else X_test
    
    print(f"-> Fitting model: {name}...")
    model.fit(X_tr, y_train)
    y_pred = model.predict(X_te)
    
    r2 = r2_score(y_test, y_pred)
    rmse = np.sqrt(mean_squared_error(y_test, y_pred))
    mae = mean_absolute_error(y_test, y_pred)
    
    results.append({
        "Model": name,
        "R2 Score": round(r2, 4),
        "RMSE": round(rmse, 2),
        "MAE": round(mae, 2)
    })

# 4. SUMMARY EVALUATION TABLE
print("\n" + "-" * 70)
print("[STEP 3] MODEL PERFORMANCE COMPARISON SUMMARY")
print("-" * 70)

summary_df = pd.DataFrame(results).sort_values(by="R2 Score", ascending=False).reset_index(drop=True)
print(summary_df.to_string(index=False))

# 5. FEATURE IMPORTANCES
print("\n" + "-" * 70)
print("[STEP 4] Top 10 Feature Importances (Random Forest)")
print("-" * 70)

rf_model = models["Random Forest"][0]
importances = pd.DataFrame({
    'Feature': feature_cols,
    'Importance': rf_model.feature_importances_
}).sort_values(by='Importance', ascending=False).reset_index(drop=True)

print(importances.head(10).to_string(index=False))

print("\n" + "=" * 70)
print("PHASE 3 COMPLETE: MULTI-MODEL BENCHMARK FINISHED")
print("=" * 70)

======================================================================
PHASE 3: MULTI-MODEL ML BENCHMARK (MSE FEATURE ENRICHED - NO DEPTH)
======================================================================
-> Successfully Engineered 'MSE_psi' and 'MSE_psi_roll5_mean' Features.

[STEP 0] Cleaned Feature Matrix Loaded: 6812 active rows.
-> Predictor Features (Excluding Raw Depth): 47

----------------------------------------------------------------------
[STEP 1] Performing Block Depth Split then Shuffling Train Set...
-> Training Set (Shuffled):   5449 samples
-> Testing Set (Sequential): 1363 samples

----------------------------------------------------------------------
[STEP 2] Training Benchmark Machine Learning Models...
-> Fitting model: Random Forest...
-> Fitting model: Gradient Boosted Trees...
-> Fitting model: Support Vector Machine (SVR)...
-> Fitting model: Nearest Neighbors (k-NN)...
-> Fitting model: Polynomial Ridge Regression...
-> Fitting model: Multilayer Perceptron (Neural Net)...

----------------------------------------------------------------------
[STEP 3] MODEL PERFORMANCE COMPARISON SUMMARY
----------------------------------------------------------------------
                             Model   R2 Score   RMSE   MAE
            Gradient Boosted Trees     0.5390   3.29  1.77
                     Random Forest     0.2782   4.11  2.24
Multilayer Perceptron (Neural Net)    -0.0589   4.98  3.74
          Nearest Neighbors (k-NN)    -0.2702   5.46  4.13
      Support Vector Machine (SVR)    -4.0890  10.92  5.30
       Polynomial Ridge Regression -2986.7756 264.68 36.75

----------------------------------------------------------------------
[STEP 4] Top 10 Feature Importances (Random Forest)
----------------------------------------------------------------------
                      Feature  Importance
  Pump Press (psi)_roll7_mean    0.573597
            Flow In (gal/min)    0.117394
                  TARGET_lag1    0.086236
                      MSE_psi    0.070967
Rotary Speed (rpm)_roll3_mean    0.050004
       Flow In (gal/min)_lag2    0.014946
 Flow In (gal/min)_roll7_mean    0.010516
  Pump Press (psi)_roll3_mean    0.007617
  Hookload (k-lbs)_roll3_mean    0.007552
             Pump Press (psi)    0.007467

======================================================================
PHASE 3 COMPLETE: MULTI-MODEL BENCHMARK FINISHED
======================================================================

Phase 5: Mechanical Specific Energy (MSE) Model Benchmark

---

Key Findings & Performance Breakthrough

By integrating Mechanical Specific Energy (MSE) into the feature set, the model achieved a major breakthrough in predictive stability without relying on Depth as an explicit feature:

• Model Stabilization: Gradient Boosted Decision Trees improved drastically from an unviable $R^2 = -1.3525$ (in the prior depth-free run) to an $R^2 = 0.5390$.
• Error Reduction: Mean Absolute Error (MAE) dropped to $1.77\text{ ft/hr}$, demonstrating strong predictive precision across unseen test intervals.
• Physical Proxy Mechanism: MSE successfully serves as a mathematical proxy for rock strength by combining mechanical work and hydraulic energy ($WOB$, $RPM$, $Torque$, and $ROP$).
• Extrapolation Ability: Incorporating surface-measured drilling physics allows tree-based algorithms to generalize into deeper, harder formations using surface operational parameters alone.

5.3 Transient Dynamics & Rolling Window Interpretation

The feature importance distribution confirms that drilling response is inherently transient rather than instantaneous:

• Hydraulic Smoothers (Pump Press (psi)_roll7_mean - (57.36%)): Surface standpipe pressure experiences time delays and dampening due to mud column compressibility and downhole pressure wave propagation. Rolling averages capture steady-state hydraulic energy better than instantaneous noisy sensor readings.
• Autoregressive Memory (TARGET_lag1 - (8.62%)): Including lagged target (ROP) introduces autoregressive temporal memory, allowing models to maintain continuity across smooth formation boundaries.

---

Core Chain of Thought (Quick Review)

The Problem: Removing Depth caused models to fail ($R^2 < 0$) because surface parameters alone couldn't sense changing rock hardness downhole.

The Solution: Calculate MSE to capture mechanical efficiency and rock strength independently of depth.

The Result: The model stabilized ($R^2 = 0.5390$, $\text{MAE} = 1.77\text{ ft/hr}$) and learned to extrapolate across new formations purely through surface drilling physics.

Phase 6: Time-Series Cross-Validation Strategy

---

Strategy & Methodological Intent

To validate model stability and prevent spatial data leakage, we implement a Time-Series Cross-Validation strategy to evaluate predictive performance on unseen sequential data.

• Mitigating Data Leakage: Standard k-fold cross-validation randomly shuffles data, which introduces data leakage in time-series drilling logs due to strong temporal and spatial autocorrelation.
• Temporal Integrity: Time-Series CV respects sequential order by training strictly on past drilling intervals to predict future formation segments.
• Generalization Assessment: This approach provides a realistic measure of how accurately the model will perform in real-time operations when encountering independent, unseen formations downhole.

6.1 Temporal Integrity & Data Leakage Prevention

In spatial logging data, consecutive measurements recorded at close depth intervals share near-identical formation characteristics and signal dynamics.

• The Leakage Risk of Random (K)-Fold CV: Randomly shuffling rows before splitting mixes neighboring spatial samples across training and testing sets. This causes models to artificially "memorize" nearby formation states, yielding overly optimistic (R^2) scores that fail in actual forward drilling operations.
• Block Depth & Sequential Validation: By enforcing a strict Block Depth Split (Training on (80%) sequential upper-hole samples, Testing on (20%) unseen deeper sequential samples), the model is forced to generalize strictly from surface mechanical and hydraulic relationships without prior exposure to the test lithology.

---

Review

The Goal: Confirm the model is not overfitted and can reliably predict on true unseen data.

The Problem: Random cross-validation cheats by mixing past and future data points together.

The Solution: Use Time-Series CV to train on early data and test on subsequent formations, accurately simulating real-time forward drilling.

import pandas as pd
import numpy as np

from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import TimeSeriesSplit, GridSearchCV
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

print("=" * 70)
print("PHASE 4: TIME-SERIES CROSS-VALIDATION & HYPERPARAMETER TUNING")
print("=" * 70)
 
# 0. LOAD AND PREPARE DATASET (MSE ENRICHED, EXCLUDING DEPTH)
file_input = "drilling_phase2_engineered.csv"
df_p4 = pd.read_csv(file_input)

depth_col = 'Depth(ft)' if 'Depth(ft)' in df_p4.columns else 'Depth'
target_col = 'ROP(1 ft)' if 'ROP(1 ft)' in df_p4.columns else 'ROP'

# Engineer MSE (Teale's Equation)
bit_diameter = 8.5
bit_area = (np.pi / 4.0) * (bit_diameter ** 2)
wob_col = [c for c in df_p4.columns if 'WOB' in c.upper() or 'WEIGHT ON BIT' in c.upper()]
rpm_col = [c for c in df_p4.columns if 'RPM' in c.upper() or 'ROTARY' in c.upper()]
torque_col = [c for c in df_p4.columns if 'TORQUE' in c.upper()]

if wob_col and rpm_col and torque_col:
    wob = df_p4[wob_col[0]]
    rpm = df_p4[rpm_col[0]]
    torque = df_p4[torque_col[0]]
    rop = df_p4[target_col].replace(0, np.nan)
    
    df_p4['MSE_psi'] = (wob / bit_area) + ((120 * np.pi * rpm * torque) / (bit_area * rop))
    df_p4['MSE_psi_roll5_mean'] = df_p4['MSE_psi'].rolling(window=5, min_periods=1).mean()

# Ensure dataset is sorted strictly by depth
df_p4 = df_p4.sort_values(depth_col).reset_index(drop=True)

# Exclude depth from predictor features
ignore_cols = [depth_col, target_col]
feature_cols = [col for col in df_p4.columns if col not in ignore_cols]

# Clean Inf/NaNs
df_p4[feature_cols] = df_p4[feature_cols].replace([np.inf, -np.inf], np.nan)
df_p4 = df_p4.dropna(subset=feature_cols + [target_col]).reset_index(drop=True)

# 1. SEQUENTIAL BLOCK SPLIT (80% Train Block / 20% Holdout Test Block)
split_idx = int(len(df_p4) * 0.80)
train_block = df_p4.iloc[:split_idx]
test_block  = df_p4.iloc[split_idx:]

X_train_seq = train_block[feature_cols]
y_train_seq = train_block[target_col]

X_test = test_block[feature_cols]
y_test = test_block[target_col]

print(f"\n[STEP 1] Data Prepared: {len(X_train_seq)} training rows, {len(X_test)} holdout test rows.")

# 2. STREAMLINED GRID SEARCH (3 Folds, n_jobs=1 for single-thread stability)
tscv = TimeSeriesSplit(n_splits=3)

param_grid = {
    'n_estimators': [100, 200],
    'learning_rate': [0.03, 0.05],
    'max_depth': [3, 4],
    'subsample': [0.8, 1.0]
}

gbm_base = GradientBoostingRegressor(random_state=42)

print("\n[STEP 2] Running GridSearchCV with TimeSeriesSplit...")
grid_search = GridSearchCV(
    estimator=gbm_base,
    param_grid=param_grid,
    cv=tscv,
    scoring='neg_mean_squared_error',
    n_jobs=1,         # Prevents Windows multi-threading CPU locks
    verbose=2         # Prints continuous step-by-step progress
)

grid_search.fit(X_train_seq, y_train_seq)

best_gbm = grid_search.best_estimator_

print("\n" + "=" * 70)
print("OPTIMAL HYPERPARAMETERS FOUND BY CROSS-VALIDATION")
print("=" * 70)
for k, v in grid_search.best_params_.items():
    print(f"   - {k}: {v}")

# 3. EVALUATE TUNED MODEL ON UNSEEN DEEP HOLDOUT TEST SET
y_pred_tuned = best_gbm.predict(X_test)

r2_tuned = r2_score(y_test, y_pred_tuned)
rmse_tuned = np.sqrt(mean_squared_error(y_test, y_pred_tuned))
mae_tuned = mean_absolute_error(y_test, y_pred_tuned)

print("\n" + "=" * 70)
print("FINAL PERFORMANCE: TUNED GRADIENT BOOSTING (HOLDOUT TEST BLOCK)")
print("=" * 70)
print(f"R2 Score : {r2_tuned:.4f}")
print(f"RMSE     : {rmse_tuned:.2f} ft/hr")
print(f"MAE      : {mae_tuned:.2f} ft/hr")
print("=" * 70)

======================================================================
PHASE 4: TIME-SERIES CROSS-VALIDATION & HYPERPARAMETER TUNING
======================================================================

[STEP 1] Data Prepared: 5449 training rows, 1363 holdout test rows.

[STEP 2] Running GridSearchCV with TimeSeriesSplit...
Fitting 3 folds for each of 16 candidates, totalling 48 fits
[CV] END learning_rate=0.03, max_depth=3, n_estimators=100, subsample=0.8; total time=   3.6s
[CV] END learning_rate=0.03, max_depth=3, n_estimators=100, subsample=0.8; total time=   7.3s
[CV] END learning_rate=0.03, max_depth=3, n_estimators=100, subsample=0.8; total time=  10.1s
[CV] END learning_rate=0.03, max_depth=3, n_estimators=100, subsample=1.0; total time=   4.4s
[CV] END learning_rate=0.03, max_depth=3, n_estimators=100, subsample=1.0; total time=   8.0s
[CV] END learning_rate=0.03, max_depth=3, n_estimators=100, subsample=1.0; total time=  12.8s
[CV] END learning_rate=0.03, max_depth=3, n_estimators=200, subsample=0.8; total time=   7.2s
[CV] END learning_rate=0.03, max_depth=3, n_estimators=200, subsample=0.8; total time=  14.5s
[CV] END learning_rate=0.03, max_depth=3, n_estimators=200, subsample=0.8; total time=  21.3s
[CV] END learning_rate=0.03, max_depth=3, n_estimators=200, subsample=1.0; total time=   8.1s
[CV] END learning_rate=0.03, max_depth=3, n_estimators=200, subsample=1.0; total time=  17.1s
[CV] END learning_rate=0.03, max_depth=3, n_estimators=200, subsample=1.0; total time=  26.5s
[CV] END learning_rate=0.03, max_depth=4, n_estimators=100, subsample=0.8; total time=   4.4s
[CV] END learning_rate=0.03, max_depth=4, n_estimators=100, subsample=0.8; total time=  10.1s
[CV] END learning_rate=0.03, max_depth=4, n_estimators=100, subsample=0.8; total time=  14.2s
[CV] END learning_rate=0.03, max_depth=4, n_estimators=100, subsample=1.0; total time=   4.9s
[CV] END learning_rate=0.03, max_depth=4, n_estimators=100, subsample=1.0; total time=  11.8s
[CV] END learning_rate=0.03, max_depth=4, n_estimators=100, subsample=1.0; total time=  16.7s
[CV] END learning_rate=0.03, max_depth=4, n_estimators=200, subsample=0.8; total time=   8.9s
[CV] END learning_rate=0.03, max_depth=4, n_estimators=200, subsample=0.8; total time=  18.9s
[CV] END learning_rate=0.03, max_depth=4, n_estimators=200, subsample=0.8; total time=  28.0s
[CV] END learning_rate=0.03, max_depth=4, n_estimators=200, subsample=1.0; total time=  10.9s
[CV] END learning_rate=0.03, max_depth=4, n_estimators=200, subsample=1.0; total time=  23.1s
[CV] END learning_rate=0.03, max_depth=4, n_estimators=200, subsample=1.0; total time=  33.9s
[CV] END learning_rate=0.05, max_depth=3, n_estimators=100, subsample=0.8; total time=   3.7s
[CV] END learning_rate=0.05, max_depth=3, n_estimators=100, subsample=0.8; total time=   6.8s
[CV] END learning_rate=0.05, max_depth=3, n_estimators=100, subsample=0.8; total time=  10.7s
[CV] END learning_rate=0.05, max_depth=3, n_estimators=100, subsample=1.0; total time=   4.0s
[CV] END learning_rate=0.05, max_depth=3, n_estimators=100, subsample=1.0; total time=   9.0s
[CV] END learning_rate=0.05, max_depth=3, n_estimators=100, subsample=1.0; total time=  14.0s
[CV] END learning_rate=0.05, max_depth=3, n_estimators=200, subsample=0.8; total time=   8.2s
[CV] END learning_rate=0.05, max_depth=3, n_estimators=200, subsample=0.8; total time=  13.6s
[CV] END learning_rate=0.05, max_depth=3, n_estimators=200, subsample=0.8; total time=  21.1s
[CV] END learning_rate=0.05, max_depth=3, n_estimators=200, subsample=1.0; total time=   8.4s
[CV] END learning_rate=0.05, max_depth=3, n_estimators=200, subsample=1.0; total time=  17.0s
[CV] END learning_rate=0.05, max_depth=3, n_estimators=200, subsample=1.0; total time=  26.5s
[CV] END learning_rate=0.05, max_depth=4, n_estimators=100, subsample=0.8; total time=   4.4s
[CV] END learning_rate=0.05, max_depth=4, n_estimators=100, subsample=0.8; total time=   9.1s
[CV] END learning_rate=0.05, max_depth=4, n_estimators=100, subsample=0.8; total time=  13.5s
[CV] END learning_rate=0.05, max_depth=4, n_estimators=100, subsample=1.0; total time=   5.9s
[CV] END learning_rate=0.05, max_depth=4, n_estimators=100, subsample=1.0; total time=  10.9s
[CV] END learning_rate=0.05, max_depth=4, n_estimators=100, subsample=1.0; total time=  17.4s
[CV] END learning_rate=0.05, max_depth=4, n_estimators=200, subsample=0.8; total time=   8.8s
[CV] END learning_rate=0.05, max_depth=4, n_estimators=200, subsample=0.8; total time=  17.8s
[CV] END learning_rate=0.05, max_depth=4, n_estimators=200, subsample=0.8; total time=  28.4s
[CV] END learning_rate=0.05, max_depth=4, n_estimators=200, subsample=1.0; total time=  11.3s
[CV] END learning_rate=0.05, max_depth=4, n_estimators=200, subsample=1.0; total time=  21.3s
[CV] END learning_rate=0.05, max_depth=4, n_estimators=200, subsample=1.0; total time=  34.6s

======================================================================
OPTIMAL HYPERPARAMETERS FOUND BY CROSS-VALIDATION
======================================================================
   - learning_rate: 0.05
   - max_depth: 4
   - n_estimators: 200
   - subsample: 1.0

======================================================================
FINAL PERFORMANCE: TUNED GRADIENT BOOSTING (HOLDOUT TEST BLOCK)
======================================================================
R2 Score : 0.7219
RMSE     : 2.55 ft/hr
MAE      : 1.53 ft/hr
======================================================================

Phase 7: Final Model Performance & Geological Validation

---

Performance Benchmarks (Tuned Gradient Boosting)

• $R^2$ Score: 0.7219 (Strong correlation and variance explained on the holdout test block)
• Root Mean Squared Error (RMSE): $2.55\text{ ft/hr}$
• Mean Absolute Error (MAE): $1.53\text{ ft/hr}$

---

Loss Function Rationale: Mean Squared Error (MSE)

We selected Mean Squared Error (MSE) as our loss function for model optimization and hyperparameter tuning:

• Quadratic Error Penalty: MSE penalizes larger prediction discrepancies quadratically.
• High-Variance Drilling Events: This forces the Gradient Boosting model to heavily penalize and minimize significant errors during abrupt downhole shifts, such as formation transitions or severe bit bounce.

---

Model Validation & Depth Track Results

To verify model performance against real-world geological conditions, we implemented two key visualizations:

1. Actual vs. Predicted Evaluation (Scatter Plot): Validates residual spread, confirming strong fit and minimal bias across performance metrics ($R^2$, RMSE, MAE).
2. Continuous Depth Track Visualization: Displays continuous actual measured ROP (Black) vs. model-predicted ROP (Orange). The high degree of overlay across depth intervals confirms that the model accurately tracks mechanical response changes caused by downhole lithology shifts.

---

Review

The Tuning Goal: Optimize Gradient Boosting hyperparameters using MSE loss to penalize large errors during sharp lithology shifts.

The Outcome: Pushed $R^2$ up to 0.7219 and reduced MAE to $1.53\text{ ft/hr}$ on unseen formations.

The Validation: Overlaying predicted ROP (Orange) on actual ROP (Black) down the depth track proves the model successfully adapts to physical formation changes in real time.

import matplotlib.pyplot as plt
import seaborn as sns

# Set clean visual styling
sns.set_theme(style="whitegrid")
fig, axes = plt.subplots(1, 2, figsize=(16, 6))

# 1. ACTUAL VS. PREDICTED PARITY SCATTER PLOT
axes[0].scatter(y_test, y_pred_tuned, alpha=0.6, color='teal', edgecolors='k', linewidth=0.5, label='Holdout Predictions')

# Add 1:1 Parity Reference Line (Ideal Match)
min_val = min(y_test.min(), y_pred_tuned.min())
max_val = max(y_test.max(), y_pred_tuned.max())
axes[0].plot([min_val, max_val], [min_val, max_val], 'r--', linewidth=2, label='1:1 Parity Line (Ideal)')

axes[0].set_title(f"Tuned GBM: Actual vs. Predicted ROP (R² = {r2_tuned:.4f})", fontsize=12, fontweight='bold')
axes[0].set_xlabel("Actual ROP (ft/hr)", fontsize=11)
axes[0].set_ylabel("Predicted ROP (ft/hr)", fontsize=11)
axes[0].legend(loc='upper left')

# 2. DOWNHOLE DEPTH TRACK (Actual ROP vs Predicted ROP along Depth)
test_depths = df_p4.iloc[split_idx:][depth_col]

axes[1].plot(y_test.values, test_depths.values, label='Actual ROP', color='black', alpha=0.7, linewidth=1.2)
axes[1].plot(y_pred_tuned, test_depths.values, label='Tuned GBM Predicted ROP', color='orangered', alpha=0.85, linewidth=1.2)

axes[1].set_title("Downhole ROP Prediction Tracking (Deep Holdout Block)", fontsize=12, fontweight='bold')
axes[1].set_xlabel("ROP (ft/hr)", fontsize=11)
axes[1].set_ylabel("Depth (ft)", fontsize=11)
axes[1].invert_yaxis()  # Inverts y-axis so deeper depth goes downwards
axes[1].legend(loc='lower left')

plt.tight_layout()
plt.show()

Model Visualizations & Formation Adaptation

We plotted actual ROP vs Predicted ROP to show the accuracy of our model and also the Depth track to visualize how well the prediction of our model (Orange) overlays the actual measured ROP (black) as depth increasesspan_2span_2.

• Key Takeaway: This shows how successfully the model adapts to physical formation changes downholespan_3span_3.

  Feature Importance & Physical Driver Analysis

To confirm that predictive power stems from genuine engineering signals rather than depth memorization, we evaluated feature gain across the tuned tree ensembles:

• Mechanical & Energy Signals: Engineered features like Mechanical Specific Energy (MSE_psi) and its moving average (MSE_psi_roll5_mean) rank among the highest drivers of ROP prediction.
• Hydraulic & Dynamic Interaction: Standpipe pressure and rotational surface variables drive the model's ability to adapt dynamically to lithology shifts downhole without relying on explicit depth coordinates.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# Extract feature importances from the tuned Gradient Boosting model
importances = best_gbm.feature_importances_

# Sort top 15 features
feature_imp_df = pd.DataFrame({
    'Feature': feature_cols,
    'Importance': importances
}).sort_values(by='Importance', ascending=False).head(15)

# Plot Feature Importances
plt.figure(figsize=(10, 5))
sns.set_theme(style="whitegrid")

ax = sns.barplot(
    data=feature_imp_df,
    x='Importance',
    y='Feature',
    palette='Blues_r',
    hue='Feature',
    legend=False
)

plt.title('Top 15 Most Important Features in Tuned ROP Model', fontsize=12, fontweight='bold')
plt.xlabel('Relative Feature Importance (Gini/Gain)', fontsize=11)
plt.ylabel('Feature Name', fontsize=11)

# Annotate exact values on each bar
for p in ax.patches:
    width = p.get_width()
    ax.annotate(f'{width:.3f}', 
                (width + 0.002, p.get_y() + p.get_height() / 2.),
                ha='left', va='center', fontsize=9, color='black', xytext=(5, 0),
                textcoords='offset points')

plt.tight_layout()
plt.show()

Phase 8: Operational Deployment & Real-Time Application

The developed (ROP) predictive pipeline provides significant practical utility during real-time drilling operations:

1. Drilling Dysfunction Detection: Deviations where measured (ROP) falls significantly below predicted (ROP) (under constant (WOB) and (RPM)) serve as an automated indicator for bit balling, severe cutter wear, or unexpected formation shifts downhole.
2. Parameter Optimization: By evaluating predicted (ROP) against varying combinations of surface (WOB) and (RPM), the model can assist directional drillers in maintaining optimal Mechanical Specific Energy (MSE) in real time.